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We survey some mathematical aspects of finite element methods for 
incompressible viscous flows, concentrating on the steady primitive variable 
formulation. We address the discretization of a weak formulation of the Navier- 
Stokes equations; we then consider the div-stability condition, whose 
satisfaction insures the stability of the approximation. Specific choices of 
finite element spaces for the velocity and pressure are then discussed. 
Finally, the connection between different weak formulations and a variety of 
boundary conditions is explored. 
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MATHEHATICAL ASPECTS OF FINITE ELEMENT 
METHODS FOR INCOMPRESSIBLE VISCOUS FLOWS 


One of the most successful and well developed mathematical theories 
concerning finite element methods is that connected with incompressible flow 
problems. The success of this theory lies not only in the accumulated elegant 
mathematical results, but also in its impact on practical computations. The 
outstanding monographs by Girault and Raviart CGR1,GR2] give a rigorous account 
of this theory and to this day remain the definitive sources. 

In this survey we examine certain mathematical aspects of finite element 
methods for the approximate solution of incompressible flow problems. Our 
principal goal is to present some of the important mathematical results which 
are relevant to practical computations. In so doing we also discuss useful 
algorithms. Due to space limitations we focus on the steady primitive variable 
formulation. Moreover, even within this narrow context, we will concentrate on 
only one of the very many different known approaches. Some other approaches are 
discussed in, e.g., CGRl,GR2,To] . 

We state at the outset that we make no attempt at being comprehensive in 
our coverage or in our attributions. To anyone who takes offense, we sincerely 
apologize. 


I - The Primitive Variable Formulation 


Let denote a bounded, possibly multiply connected, domain in 



d=2 or 
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3, and let f denote Its boundary. As a prototype for incompressible flow 
problems we consider the Navier-Stokes equations 

(1.1) u • grad u + grad p = v Au + f in fl 
together with the incompressibility constraint 

(1.2) div u = 0 in fl 
and the boundary condition 

(1.3) u = 0 on r 

where u is the velocity field, p the pressure, f the given body force, and u 
the given constant kinematic viscosity. In (1.1) the constant density has been 
absorbed into the pressure. Whenever u and p represent non-dlmensionalized 
variables, then v is the inverse of the Reynolds number Re. 

Following our detailed discussion of the approximation of solutions of 
(1.1)-(1.3) by finite element methods, we will consider other incompressible 
flow formulations, especially as they concern boundary conditions other than 

(1.3) . 

I.l - function spaced, noj'mi^ and f.oJ^mi> 

In order to introduce a Galerkin type weak formulation through which a 
finite element approximation is determined, we first need to define some 
function spaces, associated norms and forms involving functions belonging to 
those spaces. Lucid and more detailed accounts concerning these spaces may be 
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found in, e.g., CBaC,GR2,0R] . 

2 

First we denote by L (Cl) the space of functions which are square integrable 
over (3 and which is equipped with the inner product and norm 

r 1/2 

<P,Q> = J PQ and = (q,q> 

respectively. We then define the constrained space 

= f q € L^((2> I f q = 0 ] . 

0 ^ 

2 

Thus Lq(Q) consists of square integrable functions with zero mean over Q. This 
space is used in connection with the pressure; such a constraint is needed 
since it is clear from (1.1)-(1,3) that the pressure can be determined only up 
to an arbitrary constant. Other constraints, e.g., fixing the pressure at a 
given point, may be used instead without effecting any appreciable change in 
the results discussed below. Next we define the Sobolev spaces 

H^(Q) = ( q e L^((2> I D®q e L^((2> for s=l,...,k ) 

where D® denotes any and all derivatives of order s. Thus (Q) consists of 
square integrable functions all of whose derivatives of order up to k are also 
square integrable. H (Q) comes equiped with the norm 


:iq3 




1/2 


where the summation extends over all possible derivatives of order k or less. 

0 2 1 
Clearly H (Q)-L (Q) , Of particular interest is the space H (Q) consisting of 
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functions with one square' Integrable derivative and the subspace 
HgCn> = ( q 6 1 q = 0 on r ) 


whose elements have one square integrable derivative over (2 and which vanish on 
the boundary f. These spaces have the associated norm 


(1.4) 


/■ 2 dq 2 

iqii^ = ( "qi'o ^ ) 

i=l i 


We note that for functions belonging to the semi-norm 


(1.5) 


. d 9q 2 ^1/2 

^ 2 ''ax 'o ) 

i=l i 


is actually a norm equivalent to (1.4) and thus, for such functions, (1.5) may 
be used instead of (1.4). 

For vector valued functions we use the spaces 


and 


H^(fl) 

= = 1 

[ V 1 v^ 6 H^(Q) for 1=1, ... ,d ) 

hJ(0) 

= = 1 

[ V 1 v^ e Hq(( 3) for 1=1, . . . ,d ] . 


For example, (Q) consists of vector valued functions each of whose components 
belongs to H^(fl). H^(fl) is equiped with the norm 


, d 2 ^i/2 

!v.^ = ( r ) 

1=1 


alternately, has the 


norm 


Ivl 


1 


( t iv,' 

1=1 




2 0 2 d 

Also, the inner product for functions belonging to L- (Q)~H (C3)=£L (Q)2 Is 
given by 


(u,v) = J u»v 
Q 


where there is no ambiguity possible resulting from using the same notation for 
both the inner product of scalar and vector valued functions. 

We now define the bilinear forms 


1 


(1.6) 

a(u, v) 

= V / gradu 

: gradv 

for all u 

, v€H ( n ) 



fl 



0 

and 






(1.7) 

b(v,q) 

= -J qdivv 

for 

all v€hJ(0) 

2 

and qeL^di), 



(2 





and the trilinear form 


(1.8) ciw,u,v) = r w-gradu*v for all u, v,weH^(»2). 

a 


In (1.6) and (1.8) we have that (gradu)^j=3uj/3x^ and 


gradu:gradv 


E 


1. J=1 


3u. 3%’ 
1 1 

3x . 3x . 
J J 


and w»gradu*v 


E 


1. j=l 
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Using the bilinear form b(*,»), we can define the subspace 

J Z = ( V 6 I b(v,q) = 0 for all q e L^(12> ) 

I 

1 

which consists of (weakly) dlm/tgenca (.anctLana, i.e., functions whose 

2 

divergence is orthogonal to all ^^((2) functions. Certainly any divergence free 
function, in the strong sense, belongs to Z. 

1.2 - J HaieMin ti^pe mak f-ormulat ion 

The most commonly used weak formulation of (l.i)-(1.3) is the following. 

! 

Given fel? <(1) , we seek ueH^(l2) and P6 Lq(< 2) such that 

(1.9) a(u,v) + c(u,u,v) + b(v,p) = (f,v) for all veH^(Q) 

I 

I (1.10) b(u,q) = 0 for all qeL^(Q). 

By virtue of (1.10) we see that the solution u belongs to Z. 

2 

1 We note that L (fl) is not the largest function space for the data f such 

I 

that the problem (1.9)- (1.10) makes sense; indeed, all that is required of the 

} data is that the right hand side of (1.9) be bounded and this is possible for 

j some functions which are not square integrable. However, for our purposes, 

■ 2 

I feL (Q) is sufficiently general. 

I 

i It can be easily verified that whenever a pair u,p satisfies (1.9)-(1.10) 

and is sufficiently smooth to allow for the appropriate integrations by parts, 

I 

then u,p is also a solution of (1.1)-(1.3). Of course, (1.9)-(1.10) admit 
solutions which are not sufficiently smooth to be solutions of (1.1)-(1.3); 
hence the terminology maA l-artmtat Lan and qenejuxi. ifed solution are applied to 


(1.9)-(1.10) and their solution, respectively. On the other hand, it is also 
clear that any solution of (1.1)-(1.3), i.e., a i^tjtonq. Mlation, satisfies 

<1.9)-(1.10). 

For the weak formulation <1.9)-(1. 10), the boundary condition (1.3) is an 
ejyAentiat one, i.e., it must be imposed on the candidate solution functions. 
Below, in section IV. 3, we will discuss the natural boundary conditions 
associated with the weak formulation <1.9)-(1.10) . 

We will not enter into details concerning the existence, uniqueness, 
continuous dependence on data and regularity of solutions of (1.9)-(1.10). Such 
results may be found in, e.g., the definitive treatise of Teman CTel. 
Furthermore, many of these results are similar to those discussed below for the 
approximate problem. 


II “ The Finite Element Problem and the Divest ability Condition 

II. 1 - TAe dii>cretz finite element problem 

Once the Galerkin formulation (1.9)-(1.10) is established, the approximate 

problem which determines the finite element solution is defined in the usual 

manner. First one chooses the approximating finite element spaces, or more 

precisely, a family of finite element spaces, V^ and for the velocity and 

pressure, respectively. Here h is a parameter which is usually related to the 

size of the grid associated with the finite element partitioning of Q. Then one 

requires that (1.9)- (1.10) hold for functions belonging to these finite 

fl Jl tl tl 

dimensional spaces, i.e., one seeks u eV^ and p eS such that 

, h h. , , h h , h h ... h. . ,, h _h 

a(u ,v ) + c(u ,u ,v ) + b(v ,p ) = (f,v ) for all v elT 


( 2 . 1 ) 
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and 

( 2 . 2 ) 


b(u^,q^) = 0 


for all q^eS^. 


If and are subspaces of the underlying infinite dimensional spaces of 


(1.9)- (1.10), i.e., if 7^cH^(fl) and S"cLq(Q), then the finite element solution 

defined by (2.1)-(2.2) is said to be c<3nf.<3jmLng.. Otherwise, i.e., if (Q) 

h 2 

and/or S ^L^({2), then the method is said to be nan-canf.ajming.. We will restrict 
our attention to e>camples of the former. 

Once one chooses specific bases for and S^, (2.1)-(2.2) are equivalent 
to a nontlnear a(. alqjzlijviic e(iao.tLoni>, Indeed, if £qj(x)}, j=l,...,J and 

{Vj^(x)l, k=l,...,K, denote bases sets for and V^, respectively, we may then 
write 


h , 2 , 


. J ^ K 

P = Y. ^ 

J=1 k=l 


for some constants a., j=l,...,J, and B., k=l,...,K. Substituting into (2.D- 

J k 


(2.2) then yields 


k=l k,m=l 


m 


(2.3) 


+ ^ b(v^,qj) 

j=l 


and 
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K 

(2.4) for 

k=l 


which constitute a nonlinear algebraic, in fact, quadratic, system of J+K 

equations for the J+K unknowns a., and 6,, k=l,...,K. Note that the 

J ^ 

the discrete continuity equation (2.2) yields the JxK rectangular linear system 

(2.4). 

II. 2 - The d Lit- ah LI it 1 / c and it Lon 

In the positive definite case, e.g., for the equations of linear 
elasticity, the mere inclusion of the finite element spaces within the 
underlying function spaces is essentially sufficient to assure that the 
approximations are well defined and are as accurate as possible for the type of 
finite elements functions being used. Here the inclusions and S^cL^(Q) 

are not by themselves sufficient to produce stable, meaningful approximations. 
We find ourselves in the realm of what are known as mixed f^Lnite element 
methods. 

There are number of conditions which the elements belonging to the finite 
element spaces should satisfy. Most of them, e.g., the boundedfiess of the 
various bilinear and trilinear forms, are easily satisfied by conforming finite 
element spaces. The one condition which presents a problem has the following 
mathematical realization: 


given any q^eS^, 


sup 


b<v^,q^) \ 
I'v^" ^ 




r'lq 


0 


(2.5) 
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where the constant y>0 may be chosen Independent of h and of the 

particular choice of q^eS^. 

This condition may be equivalently expressed in the form: 

given any q^eS^ there exists a non-zero such that 

(2.6) b(v^,q^> i yllq^’il^llv^ll^ 

where the constant y>0 may be chosen independent of h and of the 

particular choice of q^eS^. 

Of course, for each q^ a different v^ may be used in order to satisfy (2.6). 

The condition (2.S), or equivalently (2.6), is variously known as the 

Zadi^fhen:i>kaffa-%3./yu^fia-2jvffi or the 3KBS or the Ln(.-dup condition, the latter 
designation following from the third equivalent form: 


there exists a y>0, independent of h, such that 


(2.7) 


. , h h. 
f b(v ,q ) 

inf sup — r r — 

Oxq^S^ '• llv^ll^I!q"l!Q 


^ y . 


We will refer to any of the equivalent statements (2.S)-(2.7) as the condition 
for div-i^tabil it\f. Note that these have nothing to do with the non-linearity of 
the Navier-Stokes equations and, in fact, the possible problems its 
satisfaction poses is shared by the linear equations of Stokes flow. 

Associated with the finite element spaces and and the bilinear form 
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b(*, •) we have the subspace 

= ( V 6 7^ I b(v^,q^) = 0 for all e ) . 


of ditfe/f^nce ^/<ee ^unct ian^. In general Z dZ, even when 7 cHq(Q) 

h 2 

and S cLq(Q), i.e., discretely solenoidal functions are not necessarly 
solenoidal. This is, of course, entirely analogous to the finite difference 
case, e.g., a function satifying a difference approximation to the 
incompressibility constraint is not in general solenoidal. A measure of the 
"angle" between the spaces Z^ and Z is given by 


( 2 . 8 ) 


9 = sup inf I!z - z^ii^ . 

z^Z^ zeZ 

llz^ll =1 
1 


In general, 0^0^!, which is easily seen by observing that for Z^eZ, 0=0, and 
that by choosing z=0, 0=1. 

fl h 

Note that because of (2.2), the approximate velocity ueZ”, i.e., u is 
discretely solenoidal. However, since in general 2^dZ, divu^/0. Loosly 
speaking, the div-stability condition (2.S) ensures, as h-*0 at least, that 
discretely solenoidal functions tend to solenoidal functions. 

II. 3 - and cancej^ning. th& app-raximate Aolat ian 

We now present some of the available mathematical results concerning the 
solution u^,p^ of the finite element problem (2.1)-(2.2). Here we a^Mum that 
the chosen finite element spaces 7^ and satisfy the div-stability condition 
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(2.5). Subsequently, we will look into the issue of verifying that condition. 
The summary presented is based on the detailed analysis found in CCR, JR, GRl, 
GR2, and GP]. 

2 h h 

First off, for any feL^dl), (2.1)-<2.2) has a solution u ,p , provided that 
the div-stability condition (2.5) holds. However, one can prove that the 
solution is unique only for "sufficiently small" data f or "sufficiently large" 
viscosity v. More precisely, let 


K 


h 


u 



■ h h h. 
a(w ,u ,v ) 

lu^l . I v^l. Iw^l . 
Ill 






For standard choices of finite element spaces < can be shown to be independent 
of h and, in fact, depends only on (2cR*^ and d. Then, one can show that (2.1)- 
(2.2) has a unique solution whenever 


r \c h 

, J 


sup 


(2 


I h, 

Iv 1, 


i 1 . 


This condition is very similar to the one which is needed to show the 
uniqueness of the solution of (1.9)-(1.10) and in fact the latter implies the 
former, i.e., whenever (1.9)-(1.10) can be shown to have a unique solution, 
then, provided the div- stability condition is satisfied, (2.1)-(2.2) also have 
a unique solution. 

When one can show that (1.9)-(1.10) has a unique solution, it can also be 
shown that the finite element solution of (2.1)-(2.2) converges to that 
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solution. In addition, something can be said about the convergence of the 
finite element solution even when (1.9)- (1.10) does not posses a unique 
solution. For details, see CGRl and GR2]. 

Error estimates can also be derived. Provided that the div-stability 
condition is satisfied, we have that 

(2.9) llu - u^ll. < C. inf llu - v^ll . + C_9 inf !lp - q^II. 

V 6lr q eS 

and 

(2.10) lip - p^!l. i C- inf llu - v^!l. + inf lip - q^^' 

0 ^ h ^ 1 ^ h 0 

V eV q eS 

where 8 is defined in (2.8) and C^, i=l,..,4, are constants independent of h. 

These estimates are optimal for the "graph norm" 'lull ^+!1 pIIq of functions 

belonging to H^(Q)xLq(( 2) in the sense that the rate of convergence of the 

finite element solution, measured in this norm, is the same as that of the best 

approximation to u and p out of 7^ and S^, respectively. 

If the solution of (1.9)-(1.10) , or more precisely, of the* linearized 

adjoint problem corresponding to (1.9)-(1.10) , is sufficiently regular, then 

2 

one can obtain an improved velocity error estimate in the L (Q)-norm, namely 

(2.11) llu - i C,,hllu - u^IL 

0 5 1 

where again is independent of h. 

We see that once the div-stability condition is satisfied, the error in the 
finite element approximation depends only on the ability to approximate in the 
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chosen finite element subspaces. In general, (2.9)-(2.10) indicate that the 
velocity and pressure errors are coupled. Furthermore, one finds that it is 
efficient to equilibrate the rates of convergence of the two terms on the right 
hand side of (2.9)-<2.10) . For this reason, one would like to use, for example, 
polynomials of one degree higher for the velocity components than those used 
for the pressure. As a final comment, we note that the constants appearing in 

(2.9)-(2.10) are in general proportional to 1/y where y is the stability 
constant appearing in (2.S). 

II. 4 - Veplf-^ing. the div-^tabit candit ian 

For particular choices of and S", it is usually not an easy matter to 
verify that the div-stability condition holds. To accomplish this task for 
families of such spaces is even more difficult. Here, we sketch three 
techniques for verifying the div-stability condition. 

a) Fortin’s method - One seemingly attractive method of showing that the 
div- stability condition holds is due to Fortin. He has shown [FI that the div- 
stability condition (2.5) is equivalent to the existence of a linear operator 
71 ^ from such that given any veH^(Q) 

b(R^v,q^) = b(v,q^> for all q^eS^ 

and 

m\w. i c;iv!i. 

1 1 

where the constant C>0 may be chosen independent of h and of the particular 
choice of veH^(Q). Thus the task of verifying the div-stability condition (2.5) 
is reduced to the task of showing the existence of the operator /I^ ; 
unfortunately, although the latter task has been accomplished in a few specific 
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settings, in general, it is also a difficult thing to do. 

b) VerfOrth's method - Verfurth CVll has developed a method for verifying 
the div-stability condition (2.5) which applies to the case of cant inuaui> 
discrete pressure spaces. 

Specifically, if S^cH^ (fl)riLQ((2), he starts out by combining the inverse 
inequality, see, e.g., CCil, 

(2.12) ^ . 

and the result 

b(v^,q^) h h h 

(2.13) sup r C_Iq 1. for all q eS 

h h .. h|, 2^1 

0/v"eV" V II Q 


to yield 


(2.14) 


sup 


h_yh 


h h. 
b(v ,q ) 

, h, 

Iv 1, 


^ p hlq^i 


for all q^S^ 


The inequality (2.13) can be shqvm to hold for many element pairs involving 
continuous discrete pressure fields; see, e.g., CBPl. Note that (2.13) has a 
similar appearance to the div-stability condition (2.5), but that it involves 
the "wrong" norms. 

Next, one combines the result, which can be found in, e.g., LGR1,GR2, L] ; 
given any q^eS^cL^(^2> , there exists a we 1^(Q> such that divw=q^ and 
Iwl^^C^Ilq^il^, with the approximation theoretic assumption concerning the space 
: for any there exists a w^eV^ such that 
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(2. 15) for k=0,l , 
to yield 

(2.16) 

Verftirth then shows that the div-stabillty condition (2.5) follows from (2.14) 
and (2.16) provided the constants C. ,...,C* are independent of h. 

Thus the main task of applying his method, once the inverse inequality 
(2.12) and the approximation theoretic result (2. IS) have been shown to hold 
for the discrete velocity space V^, is to show that (2.13) is valid. 

c) The Boland-Nicolaides method - A more useful method, in the sense of 
having wide applicability and relative ease of use, has been developed by 
Boland and Nicolaides CBNll. One difficulty with verifying the div-stability 
condition (2.S) is its g.tabat nature; Boland and Nicolaides have shown how to 
local if e the difficult part of the verification process. 

Specifically, consider a subdivision of Q into disjoint macj^o-clermntiy 
r=l,...,R, each of which consists of one or a few elements in the the finite 
element triangulation associated with and S^. The number of elements within 
a macro- element is independent of h, i.e., as we refine the mesh the macro- 
elements are also refined so that they always contain the same number of 
elements. Let denote the boundary of the macro-element Q^. 

Now, first suppose that the div-stability condition holds for the pair 

h ^ 

and S lacallj^ over a macro-element, i.e., there exists a constant y>0, 

independent of h and of the particular choice of macro-element, such that 


. , h h. 
b(v ,q > 

sup . — ;; — ^ 


Cghlq 


for all q^eS^ with !lq^il^=l. 


Oi<V 6l 
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(2.17) 


f b(v^,q^) ^ h K h 

gyp I __:i_ J ^ yl|q"il^ for all q es" 

11 V IL 

r 1 




where 


^ = ( veV ^ I I v=0 on ) and sjl = ( qeS^Ijj ) 

r* 


Since v|| and s|| have fixed small dimension, independent of h, (2.17) may often 
be. verified by a direct computation. 

Second, suppose that the div-stability condition holds for the 


spaces ^ and S" where 


sh 


(2.18) 


I V^cV^ 

1 f 2 

I _ ^0^^^ piecewise constant functions with respect "j 

to the macro-elements fl , r=l,...,R j ’ 


i.e., suppose that there exists a constant f>0, independent of h, such that 


(2.19) 


sup 
h Tilt (■ 


, , h h , 
b(v ,q ) 




j i for all q^eS^ . 


Sumarizing the Boland-Nicolaides method, suppose we know that the pair 
V^,S^ is larcallif cUv-^taMte with constant y independent of h, i.e., in the 
sense of (2.17). Further, suppose that the campopii>an spaces ^,S^, which 
satisfy (2.i8), are g.lobai.ti^ ditr-^toMle with constant f independent of h, i.e., 
in the sense of (2.19). Then the spaces V^,S^ are g-labailyt div-^table with a 
constant y independent of h. Thus, through the use of the comparison spaces the 
div-stability of the pair 7^,S^ need only be checked locally, i.e., over a 


macro-element . 
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This method has been succesfully used, e.g., in CBC,BN1,BN3,GR23 , to show 
the div-stability of a variety of well known elements and some novel ones as 
well, both in two and three dimensions. 

I 

II. 5 - of. unotaJbte opaceo including, the hit ineojt-conotant pair- 

There are different ways in which arbitrarily chosen finite element spaces 

i 

may fail to satisfy the dlv-stabllity condition. Here we discuss some of these 

and then give specific examples, focusing on the much studied and much 

misunderstood bilinear velocity-constant pressure pair. 

The most catastrophic type of failure is for (2.2), or equivalently (2.4), 

to imply that u^=0, l.e., the only discretely solenoldal field belonging to 7^ 

is the zero vector. The approximate solution is useless since, of course, u^=0 

cannot be a good approximate solution of the Navier-Stokes equations. This type 

of situation can usually be detected by a counting argument, i.e., the discrete 

I divergence matrix b(v.,q.), j=l,...,J and k=l,...,K, appearing in (2.4) has 

^ J 

j more independent rows than columns. 

i 

I Less catastrophic is the situation wherein for one or a few, but not all, 

q^eS^ we have that b(v^,q^)=0 for all v^eV^ so that y=0 in (2.5). This kind of 
i failure of the div-stability condition is usually easy to detec?t since it 

f 

results, in practice, in the discrete divergence matrix being rank deficient. 

' Furthermore, if these type of pressure modes q^ are the sole reason for the 

invalidity of (2.5), one may often still obtain, through a filtering process, 

■ useful approximations. 

! 

^ A more subtle failure of the div-stability condition is the case where for 

' at least some q^eS^ 

i 
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(2.20) C.hllq^Ik i sup 

for some constants and independent of h. In this case r=0(h) where r is 
the constant appearing in <2.S). In practice this may result in a loss of 
accuracy, especially for the pressure approximations. Such instabilities are 
harder to detect because, of course, computations are usually carried out using 
a finite value of h. In particular no problems such as those caused by rank 
deficient approximations to the continuity equation are encountered. This type 
of situation points out the dangers of calculating on only one grid and of not 
at least performing serious mesh refinement studies. It also points out the 
usefulness of rigorous results concerning the stability, or lack thereof, of 
finite element spaces. 

“ ^-Unstable linear-constant pair - An example of the first and most 
catastrophic instability is the following seemingly natural choice for the 
velocity and pressure finite element spaces. Let fl be a square which is 
triangulated as in the figure below. For the velocity approximations we choose 
piecewise linear functions with respect to the 
given triangulation which are continuous over Q 
and which vanish on F. For the discrete pressures 
we choose piecewise constant functions with 
respect to the same triangulation and having zero 
mean over (2. Clearly V^cH^(fl) and S^cL^(Q). For 
this choice the only discrete velocity field 
satisfying the discrete incompressibility 
constraint (2.2) is u^=0, l.e., the only 




discretely solenoidal velocity field is the zero 
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vector! One easily sees that if there are N cells to side, that the number of 
equations in (2.4) is J=dimCS^)= 2N^-1 which is greater than the number of 
columns K=dim(V^)=2(N-l)^. 

In the above example we see that the discrete incompressibilty condition 
(2.2) imposes too many constraints relative to the available velocity degrees 
of freedom. In fact, dim(S*^)>dim(V^) . In order to remedy the situation one 
must, at least, increase the dimension of relative to that of S^. 

b) The bilinear-constant element pair - We next consider the bilinear 
velocity- constant pressure pair which is often refered to as the Q^~^q element 
pair. Again consider the case of Q being a square and consider the 
"triangulation" of the figure below. We now choose to consist of piecewise 
bilinear functions with respect to this triangulation which are continuous over 
Q and which vanish on F. For we choose piecewise constant functions over 
the same triangulation and which have zero mean over Q. Once again the 
inclusions and S*’cLq(Q) hold. The simple counting argument used for 

the first example does not yield any definitive information since dira(F^)= 

2 h 2 

2(N-1) , the same as before, while now dim(S )=N -1. 

It is well known, e.g., see IF, BH, SGLGE, 

JP, GNP] , that this bilinear-constant element | I * I 

pair exhibits the disastrous "checkerboard" mode, 
i.e., for the particular discrete pressure field 
q^eS^ which is +1 in the "red squares" and -1 in ' 

the "black squares" we have that b(v^,q^)=0 for 

all v^eV^. This is an example of the second type 
of instability discussed above. The single "bad" 


pressure mode can be easily filtered out, and 
therefore some have suggested that once this mode 
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is taken care of, the bilinear-constant element pair can be safely used. 

However, this is not the whole story for the bilinear-constant element 

pair. Boland and Nicolaides CBN23 have shown that there exist other pressure 

modes for which (2.20) is satisfied. The left hand inequality of (2.4) was 

previously known CJPl, at least in the different context of penalty methods. Of 

course, the left inequality does not imply the right, and certainly doesn’t 

imply that for those modes the stability constant y=0(h>. However, Boland and 

Nicolaides have shown that this is indeed the case. Moreover, they have shown 

[BN3] that there exist data f for which the pressure approximations do not 

converge and that it is also possible to set up problems for which the velocity 

approximations do not converge as well. At the least, since the constants in 

-1 

the error estimates (2.9)-(2.11) are proportional to y , there will likely be 
a loss of accuracy due to these pressure modes. Their conclusions are worth 
noting, especially in view of the fact that the bilinear- constant element 
pair, with the checkerboard mode filtered out, has been used on numerous 
ocasslons in "practical" computations. 


Ill - Finite Element Spaces for the Primitive Variable Formalation 

In this section we discuss pressure and velocity finite element spaces 
which have been rigorously sho\m to satisfy the div-stability condition. There 
are many such pairs known, especially for two dimensional problems; therefore 
we will restrict our attention to pairs which have proven to be of the most 
practical utility. 

Throughout, denotes the space of polynomials of degree less than 

equal to k with respect to the set SDcR^ and CP^(3))1^ denotes the space of d- 
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vector valued functions each of whose components belong to P^(3)). Analogous 
definitions hold for and in the case of functions which are 

polynomials of degree less than or equal to k in each of the coordinate 
directions, e.g., Q^(^> denotes piecewise bilinear functions with respect to 
the set 2). Likewise we define the spaces C^(3D) and of k-times 

continuously differentiable functions with respect to the set 3). 

2 

For the most part, the results below hold for polygonal domains in R and 

3 

polyhedral domains in R . Through the use of, e;g., isoparametric elements, 
they will also hold for domains with curved bouhdaries provided the latter 
satisfy the usual smoothness criteria. Furthermore, we assume that all 
subdivisions of Q into finite elements which are employed below satisfy the 
standard . conditions. For details concerning these issues, one may consult, 
e.g., CCi], 

Ill.i - ieceu/i^ linecLr and bltlmar ueioclt^ ^letd^ 

We begin with some examples involving piecewise linear or bilinear velocity 
fields with respect to a subdivision of ^2 into triangles or rectangles, 

respectively. In all cases the discrete velocity fields are continuous over Q. 

' « 

In combination with these type of velocity finite element spaces we will 

consider both discontinuous piecewise constant and continuous, over Q, 
piecewise linear pressure fields. Every element pair listed satisfies the div- 
stability condition (2.5). Moreover, provided the solution u,p of (1.9>-(1.10) 
satisfies and psH^ (Q)f)L^(Q) , the following error estimates for 

the discrete solution u^,p^ of (2.i)-(2.2) hold uniformily in h: 

, - u^:i . = 0(h) 

I 

j "u - u^l'^ - O(h^) 

'■P - = 0(h) 


( 3 . 1 ) 
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Thus, these elements yield first order accurate pressure approximations and 
second order accurate velocity approximations. 

a) Piecewise constant pressures I - For the linear-constant element pair 
mentioned in section II. 5 the discrete continuity equation overconstrained the 
approximate velocity field. However, by employing different grids for the 
pressure and velocity fields, the linear-constant element pair may be made 
stable. For example, consider a given triangulation of a polygonal domain >2 
into triangles. Then divide each triangle in into four triangles by joining 
the midsides, thus defining a refined triangulation example is 
provided in the figure below. 




A pressure triangle in 
Now define 


The four associated velocity 
triangles in 


r = ( q I qeP^(A) ,AeT^ ; J q = 0 ) 

(3.2) I 

I 7^ = ( V I veCP^(d)]^ , ; ve[C^(( 2)J^ ; v=0 on F ) 


so that the pressure is sought among piecewise constants with respect to the 
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trlangulation while the velocity is sought among continuous piecewise linear 
fields with respect to the finer triangulation finite element 
spaces defined by (3.2) are known to satisfy the div-stability condition (2.S) 
and thus yield optimally accurate approximations satisfying (3.1) 

b) Piecewise constant pressures II - For the unstable lineai — constant 
element pair of section II. S there was one velocity element for each pressure 
element; for the stable linear-constant element pair (3.2) there are four 
velocity triangles for each pressure triangle. Stable linear-constant element 
pairs may be defined wherein the ratio of discrete pressures to velocities is 
not so high. For example, let the velocity space be as in (3.2); now define 
the pressure space through the following choice of basis. For each triangle 
of we define three basis functions, namely piecewise constants which are 
unity in the shaded areas in figure below and zero in the unshaded areas. Of 
course, outside the particular triangle of the basis functions vanish as 
well. This pressure space consists of three out of the 



four possible piecewise constants associated with the four triangles in ^^^^2 
contained within a single triangle in T^. Moreover, there are essentially three 
times as many pressure degrees of freedom for this choice of as there are 
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for the choice made in (3.2). However, this element pair is also stable, i.e., 
satisfies the div-stability condition (2.5) and the error estimates (3.1). 

c) Piecewise linear pressures - One may also couple a piecewise linear 
velocity element with a piecewise linear pressure element and still satisfy the 
div-stability condition (2.5) and the estimates (3.1). Such a pair was 
introduced in EBP], analyzed there and in CV13, and is given by 

r = ( q I qeP (A) , AeT ; qeC®(Q) ; J q = 0 ) 

(3.3) Q 

I IT = as in (3.2). 

Due to the coupling between the pressure and velocity errors one cannot take 

advantage of the better approximating ability of the linear pressure space. 

Thus, insofar as the rates of convergence, this lineai — linear element pair is 

no better than the stable linear-constant element pairs. However, in practical 

calculations we have found this to be the best element combination involving 

linear velocity fields, better in the sense of giving more accuracy for useful 

values of h. Furthermore, this linear-linear element pair usually results in 

fewer unknowns, for the same grid, than do the linear-constant pairs. For 

example, suppose the pressure triangulation is given by the first figure of 

section II. 4 with N intervals on each side. Thus there are 2 n" triangles in 

2 

and the element pair (3.2) has 2N -1 pressure unknowns; on the hand, the number 

2 

of nodes in this triangulation is only CN+1) and thus the piecewise linear 

2 

pressure space of (3.3) has only CN+1) -1 degrees of freedom. Both element 

2 

pairs have 2(2N-1) velocity unknowns so that the linear-linear element pair 

2 

(3.3) has roughly N less degrees of freedom, for the same grid, as does the 
linear-constant element pair (3.2). 

d) Piecewise bilinear velocity fields - Entirely analogous to the 
triangular elements described above, we have the following elements involving 



-27- 


bllinear velocity fields with respect to rectangular elements. More general 
quadrilateral elements may be found from these through, e.g., Isoparametric 
mappings. 

We start with a subdividi vision of Q into rectangles, or more generally 
quadrilaterals. Subsequently we divide each rectangle into four smaller 
rectangles by joining the midsides, thus creating another subdivision ^ 
into rectangles. See the figure below. In all three velocity-pressure element, 
pairs 


■» 


A pressure rectangle in The four associated velocity 

rectangles in 

about to be described we choose the approximating velocity space to consist of 
piecwise bilinear functions with respect to the subdivision 
continuous over fl and which vanish on f, i.e., 

I 

(3.4) 7^ = ( V I ve[Q^(D)]^ , ; v=0 on T ) . 

For the first pressure space we choose piecewise constants with respect to 
the larger quadrilaterals of the subdivision and which have zero mean over 
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(2, i.e., 

= ( q I qeQQ(0> , ; Jq=0 ) . 

a 

As indidicated in the figure below, for the second pressure space we choose 
three 





1 


1 


out of the four possible piecewise constants associated with the rectangles 

belonging to which have zero mean over (2. Finally, the third pressure 

space consists of piecewise bilinear functions with respect to the subdivision 

2. which are continuous over (2 and have zero mean over fl, i.e., 
n 


i 


(3.5) 


= ( q I qeQ^fO) , Ue2^ ; qeC^((2) ; Jq=0 ) 

Q 


The three velocity-pressure elements Just described satisfy the div- 
stability condition (2.5) and the error estimates (3.1). Similar to the case 
for triangles and for the same reasons, the prefered element pair involving 
bilinear velocities is (3.4) coupled with (3.5), i.e., the bilinear velocity- 
bilinear pressure pair. 



-29- 


III. 2 - 7he 7a\ftaj^-Kaod element paiJ^ 

We next turn to quadratic and biquadratic approximate velocity fields. 
Suppose we have a triangulation of Q. Then, the Taylor-Hood element pair 
ITH] is defined by 


(3.6) 


J = ( V I , AeT^ ; velC(Q)l^ ; v=0 on T ) 

1 = ( q I qeF^(A) , ; qeC(Q) ; |q=0 ) . 

Q 


Note that we are now basing and on the same grid but on different degree 
polynomials, in contrast to (3.3), which uses the same degree polynomials but 
different grids. The element pair (3.6) satisfies the div- stability condition 
(2.5). Furthermore, if the solution (u,p) of (1.9)-(1.10) has the indicated 
smoothness, then the following error estimates hold 
uniformily in h: 


, liu - u^ll^ = 0(h"^^> , uerf"(fl)nHi(W) 

1 f 

(3.7) ^ Ilu - u^IIq = 0(h™> I whenever | and L m=2 or 3. 

MIp - p^IIq = 0(h™> ^ ^ peH®‘^(n)nLp(fl) ^ 

These results have been established by many authors, including CBP,V1,BN1]. We 
see from (3.7) that if ueH^(fl)ni^(Q) and peH^((2>r)L^(Q) then, in L^-norms, we 
have third order accurate velocity approximations and second order accurate 
pressure approximations. This is an improvement over any of the elements 
involving linear velocities. 

One should note that the* number of degrees of freedom, both of velocity and 
pressure type, associated with the use of (3.6) is identical to that associated 
with the use of (3.3), the most efficient linear velocity element. In fact, the 
structure of the discrete system resulting from a Taylor-Hood discretization is 
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in every way identical to that resulting from the use of (3.3). Therefore, the 
solution times for the Taylor-Hood and the linear-linear discrete systems are 
roughly the same if one uses the same pressure triangulation in both cases. Of 
course, the Taylor-Hood element pair will yield better accuracy than the linear- 
linear pair, provided the exact solution is sufficiently smooth. 

On the other hand, on the same grid, the assembly costs of Tayloi Hood will 
in general be higher since one needs to use higher order quadrature rules to 
integrate the higher degree polynomial integrands resulting from the Tayloi 
Hood element pair. For many solvers, the assembly time is overwhelmed by the 
solution time; therefore the increased assembly cost associated with (3.6) is 
not a serious drawback. Of course, this is further mitigated by the fact that 
for the same -accuracy, one may use a coarser grid for (3.6) than for (3.3). 

Summarizing, provided the exact solution is sufficiently smooth, the Taylor- 
Hood element pair, when compared to any of the linear velocity elements, yields 
better accuracy for essentially the same work, or alternately, will yield a 
desired level of accuracy for less cost. 

For rectangles or quadrilaterals we have the analogous pair 

f = f V 1 ve[Q„(D)]^ , De^. ; ve[C(Q)]^ ; v=0 on F ] 

(3.8) h 

I S = ( q I qeQj(D) , ; qeC(Q) ; Jq=0 ) 

a 

where 2^ denotes a subdivision of (3 into rectangles. This element pair 
satisfies the div-stability condition (2.5) and the error estimates (3.7). 

One may well ask if further efficiencies may be gained by using higher 
order elements, e.g., cubic velocities coupled with quadratic pressures. Here 
one needs to consider the trade-off between the increased accuracy of higher 
order elements and the increased complexity of those elements. As in other 
settings, e.g., structural mechanics, one generally finds that the optimum 



-31- 


seems to be achieved by quadratic elements. Furthermore, it is questionable 
that in general settings the exact solution of the Navier-Stokes equations is 
sufficiently smooth to enable the potential better accuracy of higher order 
elements. In our overall experience, we have found the best choice of velocity- 
pressure elements to be the Taylor-Hood element pair (3.6), or its 
quadrilateral counterpart (3.8). 


III. 3 - T life nee elements 

Ideally, one would like to choose the finite element spaces and so 
that the functions belonging to are at least discretely divergence free. 
Certainly if the functions belonging to are divergence free then they are 
discretely divergence free as well, i.e. , divv^=0 for all implies that 

V^=Z^. Such a case effects a great simplification since the velocity and 
pressure uncouple. Indeed, we need only solve 


a(u^,v^) + c(u^,u^,v^) = (f,v^) for all v\f^ 


for the discrete velocity field since in this case the term b(v^,q^> in 
(2.1) vanishes for any v^e7^=Z^. Also, since Z^sZ, note that in t;he velocity 
estimate (2.9), 9-0 so that the velocity error depends only on the ability to 
approximate in 

Unfortunately, although there are known some finite element pairs such that 
the functions in are at least locally divergence free, these have proven to 
be impractical, and we will not consider them here. We do mention that one 
obvious method of generating divergence free discrete vector fields is to take 
the curl of a piecewise polynomial field, i.e., of a piecewise polynomial 
streamfunction. One problem with this approach is that if one wants a 
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conforming velocity field, i.e., V^cH^({3), then the discrete streamf unction 

2 

field must be chosen to be continuously differentiable over Q. In R this, of 
course, necessitates the use of at least quintic streamfunctions over 
triangles, or cubic polynomials over macro-elements, e.g., the Clough-Toucher 
element. Non-conforming velocity fields can also be generated in this manner. 
See CCa,CN,GRl, and GR2] for details. 

III. 4 - 7hj^ee dirmnAcanat etementiy 

Compared to the two dimensional setting, there are known much fewer stable 
element pairs for three dimensional problems. However, there is great interest 
in this subject and therefore there has been substantial recent progress. Here 
we mention a few of the known stable three dimensional elements. 

In the first place, the three dimensional analogue of the Taylor-Hood 
element is known to be stable in 3-D; this may be shown by the methods of 
Verfhrth or Boland-Nicolaides. Specifically, we subdivide fl into tetrahedrons 
and use continuous piecewise quadratic polynomials for the velocity and 
continuous piecewise linear polynomials for the pressure. The accuracy of this 
combination is the same as in the two dimensional case. 

Next we consider linear-constant elements. Again, subdivi<de Q into 
tetrahedrons. For the pressure space we choose piecewise constants with 
respect to this initial subdivision. Now we subdivide each tetrahedron into 12 


smaller tetrahedrons 

by first 

joining the 

centroid of 

the faces 

to 

the 

vertices, and then 

joining the 

centroid of 

the large 

tetrahedron 

to 

the 


vertices and the centroids of the faces. For the velocity space we choose 
continuous piecewise linear polynomials with respect to the smaller 
tetrahedrons. 

Another stable linear-constant element pair is defined by first subdividing 
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Q into rectangular prisms, or more generally, into distortions of such prisms. 
For the pressure space we choose piecewise constants over the rectangular 
subregions. We subdivide each rectangular prism into 24 tetrahedrons by first 
drawing the two diagonals of each face, then Joining the centroid of the prism 
to the vertices and to the six intersection points of the face diagonals. 

Both these linear-constant element pairs are known to be stable and yield 
the same accuracy results as those for the two dimensional linear-constant 
pairs. See CBoCl for details. 


IV - Alternate Weak Forms and Bowidary Conditions 

In this section we examine some variants of the weak formulation (1.9)- 
(1.10), mostly from the viewpoint of how different boundary conditions may be 
incorporated into a finite element method using primitive variables. We again 
emphasize that there are many radically different weak formulations involving u 
and p which we will not be able to consider; we are restricting ourselves to 
variants of the most commonly used weak formulation. 

Before considering boundary conditions, we briefly consider an alternate 
formulation of the convection term in (1.9). 

IV. 1 - An atteynate ^arrmtat ion of. tha canvect ion term 

For the purpose of simplifying the analysis of the approximate solution it 
can be useful to introduce a slightly different weak formulation wherein the 
trilinear form c(«,*,*) appearing in (1.9) is replaced by the skew-symmetrized 
form introduced by Teman 



-34- 


(4.1) Z(w,u,v> = c(w,u,v) - c(w,v,u) ). 

One may easily verify that 2(u,u, v)=c(u,u, v) whenever divu=0 in (3 and u*n=0 on 
r, where n denotes the outward normal to T. Therefore, due to (1.2)- (1.3), it 
seems irrelevant whether one uses (1.8) or (4.1) in a weak formulation of the 
Naviei — Stokes equations. From an analysis point of view, the advantage of 

(4.1) is that 2(w,u, v)=-2(w, v,u) for any u,v,weH^(fl) while the analogous result 
for (1.8) holds only when divw=0 in (3 and one of u=0, v=0 or w»n=0 on f. 

We emphazise that, insofar as the accuracy of the approximations is 
concerned, it makes no difference whether one uses (1.8) or (4.1); we merely 
point out that many of the results concerning finite element approximations of 
solutions of (1.1)- (1.3) were first obtained through the use of (4.1). On the 
other hand, any implementation of (4.1) will result in more computational work 
than the analogous Implementation of (1.8). 

IV. 2 - SnhomogeneauA wlac it^ baandafl\f candit ian/> 

There are many different ways to treat inhomogeneous velocity boundary 

conditions. In practice, the overwhelming choice is to use the boundary 

2 

interpolant. We describe this method for polygonal domains (3cR ; entirely 
analogous ideas may be used in three dimensions and for domains with curved 
sides, the latter through the aid of, e.g.. Isoparametric elements. 

Consider the boundary condition 

(4.2) u=g on r 
and the set 

(4.3) ^g ~ ^ U6H^((2) 1 u satisfies (4.2) ] . 
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Note that Vq=H^(Q). The weak formulation which we will discretize is as 

2 

follows: seek and peL^((2) such that (1.9) and (1.10) hold. Note that the 

test function v still belongs to i.e., v=0 on f. 


V^cH^(Q) and S^cL^(Q). We denote by F^lp the restriction of to the boundary 


In order to pose our discrete problem we choose finite element spaces 

r, i.e., consist of functions defined on P and which can agree with the 

boundary values of at least one function belonging to V^. The finite element 
functions belonging to being, for example, piecewise polynomials, cannot in 
general satisfy the boundary condition (4.2); certainly, in general gs^V^lp. 
Therefore we choose an approximation to g, which we denote by g^, belonging to 
V^Ip. The most common choice for g^, and the one we consider here, is the 
interpolant of g in ^Ip- 

This choice is trivial to implement, which at least partially accounts for 
its popularity. For example, suppose is a Lagrange finite element space, 
i.e., one whose degrees of freedom are exclusively function values at points. 
Let k=l,...,F denote the usual finite element basis for V^. Let the first 

K of these basis functions be associated with interior nodes Xj^ so that for 
k=l,...,K, Vj^=0 for xeP. The remaining basis function k=K+l,...,K, are 

associated with nodes lying on P. In practical implementations there are 
more efficient node numbering schemes than the one we are using; however, the 
latter simplifies the explanations being attempted here. 

Choosing g^ to be the boundary Interpolant of g is then equivalent to 
writing 


K K 

(4.4) u (x) = gj^Vj^(x) + Y ' 

k=l k=K+l 

In (4.4) 6,, k=l,...,K, are the unknown coefficients to be determined; the 

K 
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coefficients of the basis functions associated with boundary nodes are simply 
set equal to g evaluated at the corresponding node. Note that (4.4) implies 
that 

h ^ 

g"(x) = Yi * 

k=K+l 


The contribution to u emanating from the second summation of (4.4) becomes 
part of the data of the discrete system of equations. 

Once an approximation is chosen, one may define the set 
= [ veV^ I v=g^ on r ) . 

s 


Note that is the finite element subspace of used in conjunction with 

the homogeneous boundary condition (1.3); also, clearly V^cH^ (0) is not a 

o 

subset of V . Now, the approximate problem may be defined as follows: seek 
S 

and p^eS^L^CQ) such that (2.1)-(2.2) hold for all and q^eS^, 

respectively. Again, the test functions v^ vanish on the boundary f. 

The whole discussion of the div-stability condition (2.5) carries over 

intact to the case of the inhomogeneous boundary (4.2); in (2.5) we still use 

the subspace 7^ of finite element velocity fields which vanish on the boundary. 

Results analagous to those of section 1 1. 3 can be derived in a fairly 

straightforward manner with the exception of some technicalities encountered 
2 

for the L (Q)-error estimate for the velocity approximation. See HGP, FGP, GR2] 
for details. 

In particular, if g^ is chosen to be the boundary interpolant of g in V^lp, 
then all the results, e.g., error estimates, concerning the finite element 


spaces discussed in section III are essentially still valid for the 
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inhomogeneous velocity boundary condition (4.2). Again, see CGP, FGP and GR21 
for details. 

I 

IV. 3 - j4ttey«nate baundojtj^ candit iani^ and f-armitat l 0 ni> of. the »ii>coai> tejm 

In this section we examine how different choices for the viscous term in 

i 

' (1.1) affect the natural boundary conditions of corresponding weak formlations. 

■ Some of this material can be found in CGLNl . 

I Due to (1.2), when u Is constant, the viscous term in (1.1) may be written 

in the various equivalent forms 


(4.S.1) vAu 


(4.5.2) div^ u( (grad u) + (grad u)^ ]j = 


(4.5.3) -wcurl(curlu) 


(4.5.4) u[ grad(divu) - curl(curlu) ) . 


i Although these different realizations are equivalent insofar as the partial 

differential equations are concerned, we shall see that each generates a 
different numerical method. 

If for some reason u is not constant or divu^O, then only (4.5.2) may be 
used. Indeed, (4.5.2) is the form of the viscous term which arises naturally 
I in the derivation of the Navier-Stokes equations from the principle of 

conservation of linear momentum and the Cauchy-Poisson constitutive equation. 
The other three forms (4.5.1), (4.5.3) and (4.5.4) are derived from (4.5.2) 

with the aid of (1.2) and the assumption that i»=constant. In (1.1) we have used 
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(4.5.1) only because this is the most popular choice in the literature; all of 
the results obtained so far hold equally well if one chooses (4.5.2) instead. 
As will be seen from the discussion below, (4.5.2) is, in general, to be 
prefered to (4.5.1). 

Denote two segments of the boundary f by and These segments may be 
empty, are not necessarily disjoint and, in fact, may be equal. Now, for fixed 
given functions g^ and g^, define the set 


V = f veH^ I v»n=g on F ; nxyxn=g on f 1 

g ®n n ®T T •' 

and the spaces 


and 


( veH^ I 


v*n=0 on r 

n 


vxn=0 on 



s = if r =r 

0 n 


S=L (fl) otherwise. 


where v*n denotes the component of v normal to the boundary F and 

n*vxn=v-(v*n)n is the projection of v onto the plane tangent to F. In the 

definition of we may use vxn=0 due to the relation vxn=nx (ipwxn) , i.e., 

2 

nxyxn=0 implies that nxv=0. In R , nxvxn=v*r where t is the unit tangent vector 
to r. 

Suppose that we wish to specify the boundary conditions 


(4.6.1) 

and 


u» n=g on r 
n n 


nxuxn=g on F . 

T T 


(4.6.2) 
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i.e., the normal velocity on and the tangential velocity on 

respectively. For all the weak formulations which we will consider involving 

any of the choices in (4.5), (4.6) will be e^^ntiat Jbcundoj^i^ candit ianit. Thus 

the trial solution functions u will satisfy (4.6), i.e., ueF , and the test 

S 

functions satisfy 

Consider the following weak formulation: for i=l,2,3 or 4, seek ue7 and 

S 

peS such that 

(4.7) a^(u,v) + b(v,p) + c(u,u,v) = (f,v) + d(v) for all 
and 

(4.8) b(u,q) = 0 for all qeS. 

Here, b(*,*) and c (♦,♦,•) remain as in (1.7) and (1.8), respectively, and f 
continues to denote the body force appearing in the momentum equation. The 
linear functional d(«) is given by 

(4.9) d(v) = J rv*n + | s*vxn 

r/r r/r 

n X 

where the functions r and s are additional data for the problem. In (4.9), for 
example, F/r^ denotes the complement of in f, i.e., xer/f^ implies that xef 

but Also, since v is an arbitrary test function, in direction vxn can be 

taken to be vectors spanning the tangent plane to F. 

The bilinear forms a^ (♦,•), i=l,...,4, depend on the choice made in (4.5) 
and, corresponding to the four choices possible in (4.5), are given by 

a^(u,v) = uj gradu:gradv 
Q 


(4.10.1) 
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(4.10.2) 32^0, v) = 1 J u( gradu + (gradu)^ ):( gradv + (gradv)^ ) 

^ (i 

(4.10.3) ag(u,V/' = uj (curlu)* (curlv) 

Q 

and 

(4.10.4) a^(u,v) = uJ (curlu)* (curlv) + (divu)(divv) . 

a 

In the customary manner, should u and p be sufficiently smooth, one can, 
through formal Integration by parts procedures, ascertain what differential 
equation problem the weak formulation (4.7)-(4.8) corresponds to. To begin 

with, we know that the boundary conditions (4.6) are satisfied since these are 
being required of the candidate trial functions u. We also find that the 

differential equations (1.1) and (1.2) are satisfied, where in (1.1) the 
viscous term is replaced according to (4.5), depending on which choice is made 
in (4.10). Finally, one finds the natural baundar\f candit iani> corresponding to 
the particular weak formulation. We will now discuss these in some detail for 
each possible choice in (4.10). * 

Corresponding to the paired choices (4.5.1) and (4.10.1) we have the 

natural boundary conditions 

(4.11.1) p - wn* gradu* n = r on and un*graduxn = s on f/r^. 

Unfortunately, these boundary conditions have no phi^^icat meaning.. Thus the 

choice (4.S.1), or equivalently (4.10.1), can only be used in conjunction with 
the boundary condition (4.6) specified on all of f, i.e., u given on f^=r^=r. 
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Next, consider the choices (4.S.1) and (4.10.1). The natural boundary 
conditions are then 


(4.11.2) 


j-p* 

1 


un* (gradu + 
wn* (gradu + 


(gradu)^) *n 
(gradu)^)xn 


-r 


-s 


on r/r 

n 

on r/r^ 


and 


Thus -r and -s are the normal and tangential stresses, respectively, on the 
boundaryi Then, for the choice (4.10.2), the possible combinations of boundary 
conditions at a point on the boundary f are as follows: we may specify the 
velocity, or we may specify the normal velocity and the tangential stress, or 
we may specify the tangential velocity and the normal stress. The latter 
combinations are useful, e.g., for free surface problems or at artificial 
outflow boundaries. Details may be found in CGLMl . 

The third choice (4.5.3), or (4.10.3), yields the natural boundary 
conditions 


(4.11.3) P = •" on r/r and a = s/u on r/r 

^ n X 

so that r and s are the pressure p and v times the vorticity (u=curlu, 

respectively, on the boundary. The possible combinations of boundary 

conditions are now: we may specify the velocity, or we may specify the normal 
velocity and the vorticity, or we may specify the tangential velocity and the 
pressure. The pressure is often used as an outflow condition; the vorticity is 
useful in exterior problems when matching to an inviscid irrotational flow 
since it is well known that the vorticity decays to its value at infinity 

faster than does the velocity. Again, details may be found in CGLNl . 

Unfortunately, although the boundary conditions associated with the use of 

(4.10.3) can be useful, in practice we cannot employ this particular 
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formulatlon of the viscous term. The reason for this is that the choice 

(4.10.3) requires the use of divergence free finite element velocity fields in 
order for the form ag(»,»> to be coercive on Z^. This condition is also needed 
to guarantee the stability of the approximations and, for the other three cases 
(4.10.1), (4.10.2) and (4.10.4), is trivially satisfied for any choice of 
conforming discrete velocity space. 

Fortunately, the boundary conditions (4.11.3) are appr^oximatel^ the natural 
boundary conditions associated with the choice (4.11.4). In fact, for (4.10.4), 
we have the natural boundary conditions 

(4.11.4) p - udivu = r on F/f and a> = s/v on f/f . 

n T 

The second of these is identical to the second of (4.11.3). If o is "small", 
and/or if we assume the incompressibility constraint holds up to portions of 
the boundary where the normal velocity is not specified, then (p-udivu) is 
essentially equal to p. Thus we recover, at least approximately, the first 
boundary condition of (4.11.3). 

In summary, when one has velocity and/or stress boundary conditions, one 
should use (4.11.2) in (4.7). and when one has velocity and/or pressure and/or 
vorticity boundary conditions the choice (4.11.4) is preferable. 

The discretization of (4.7)-(4.8) follows the usual procedures once one 
chooses the finite element spaces for the velocity and the pressure 
approximations. The natural boundary conditions are automatically acounted for 
by the inclusion of the linear functional d(*> in (4.7). The essential boundary 
conditions on the components of the velocity can be enforced in a manner 
analogous to that described in section IV. 2 for the case where the complete 
velocity is specified on the whole boundary. All material relating to the div- 
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stability condition <2. 5) is essentially still valid, and thus, insofar as that 
condition is concerned, the particular choices of finite elements discussed in 
section III may still be used. 

In actuality, there are very few rigorous error estimates available for 

boundary conditions other than the velocity. For polygonal or polyhedral 

domains (2, the error estimates of section II. 3 are still valid. However, for 

domains with carved baandarle^, using the type of weak formulations discussed 

here may result in a loss of accuracy. For example, for (4.10.2) with normal 

velocity and tangential stress boundary conditions, it was shown by Verfurth 

CV21 that there is a loss of accuracy due to a Babuska type paradox, i.e., the 

2 

limit of solutions of problems posed on polygonal approximations to (2ci? is not 
the solution of the problem posed on (2. Verfdrth [V31 has also shown how 

through the use of additional Lagrange multipliers on the boundary, a different 
weak formulation yields optimal accuracy. 
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